rm(list=ls())
library(clusterProfiler)
library(ggplot2)
library(enrichplot)
library(readxl)
library(ggpubr)
library(gridExtra)
library(pheatmap)
library(DESeq2)
library(RColorBrewer)
#library(biomaRt)

library(RColorBrewer)
color = colorRampPalette(rev(brewer.pal(n = 7, name =
                                          "RdYlBu")))(100)

path_barplot = function(all_go, top = 1000000) {
  path_ratio = do.call(rbind, sapply(all_go$GeneRatio, strsplit, '/'))
  path_ratio = apply(path_ratio, 2, as.numeric)
  
  bk_ratio = do.call(rbind, sapply(all_go$BgRatio, strsplit, '/'))
  bk_ratio = apply(bk_ratio, 2, as.numeric)
  
  plot_height = path_ratio[, 1] * bk_ratio[, 2] / path_ratio[, 2] / bk_ratio[, 1]
  
  cols = -log10(all_go$pvalue)
  
  df = data.frame(pathways = all_go$Description,
                  fc = plot_height,
                  'logp' = cols)
  df = df[order(df$logp, decreasing = T), ]
  if (top < length(cols))
    df = df[1:top, ]
  
  df = df[order(df$logp, decreasing = F), ]
  df$pathways = factor(df$pathways, levels = df$pathways)
  temp = df$logp
  temp1 = quantile(df$logp)
  temp[temp > temp1[3] * 5] = max(temp[temp <= temp1[3] * 5])
  temp[temp < temp1[3] / 5] = min(temp[temp >= temp1[3] / 5])
  df$logp = temp
  name_length = max(apply(df, 2, nchar)[, 1])
  ggplot(data = df, aes(x = pathways, y = fc,
                        fill = logp)) +
    geom_bar(stat = 'identity',
             width = 0.9,
             color = 'black') + coord_flip() +
    scale_fill_distiller(palette = "RdYlBu") +
    theme(
      text = element_text(size = 26),
      axis.text = element_text(size = 15 * 50 / name_length),
      legend.text = element_text(size = 15) ,
      panel.background = element_blank(),
      axis.line = element_line(colour = "black"),
    )
}




filename='/Volumes/hua_mac/research/aifantis/xufeng/manuscriipt/data/DESeq2-rev-rm_par1mr1_1-7samples-2groups_dge.MR-vs-Par.csv'
data=read.csv(filename,header=T,stringsAsFactors = F)
lfc=0.2
pval=0.01
data=data[!is.na(data$log2FoldChange) & abs(data$log2FoldChange)>lfc & !is.na(data$padj) & data$padj<pval,]

genes = data[, 1]
genome = 'org.Hs.eg.db'
eg <-
  bitr(genes,
       fromType = "SYMBOL",
       toType = "ENTREZID",
       OrgDb = genome)

gene <- eg$ENTREZID

all_go <-
  enrichGO(
    gene = gene,
    OrgDb = genome,
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 1.1,
    qvalueCutoff = 1.1,
    readable = TRUE
  )
all_go=all_go[all_go$pvalue<0.05,]



curheight = 900

barname = './figures/sFigure3G.pdf'

path_barplot(all_go, top = 20)
ggsave(barname, width = 10, height = curheight / 100)



